﻿using MatlabWork.Entity;
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace TestingPlatForm.MissionContrl.Models
{
    public class GeneratSpt
    {
        public int y = 2024;
        public int m = 9;
        public int d = 3;
        public int h = 16;
        public int mm = 1;
        public int ss = 0;

        private List<SateliteData> _sptList;

        public List<SateliteData> sptList
        {
            get { return _sptList; }
            set { _sptList = value; }
        }

        public List<double> PotoSig(List<double> XX)
        { 
            List<double> ZZ = new List<double>();
            double PI = Math.PI;
            double GE = 398600.140 * Math.Pow(10, 9);
            double Ae = 6378140;
            double D10 = 1.0;
            double D20 = 2.0;
            double X = XX[0] / Ae;
            double Y = XX[1] / Ae;
            double Z = XX[2] / Ae;
            double XD = XX[3] / Math.Sqrt(GE/Ae);
            double YD = XX[4] / Math.Sqrt(GE/Ae);
            double ZD = XX[5] / Math.Sqrt(GE/Ae);
            double R = Math.Sqrt(X*X +Y*Y +Z*Z);
            double A = D10 / (D20 / R - XD * XD - YD * YD - ZD * ZD);
            double EC = D10 - R / A;
            double ES = (X * XD + Y * YD + Z * ZD) / Math.Sqrt(A);
            double E = Math.Sqrt(EC * EC + ES * ES);
            double EE = Math.Atan2(ES, EC);
            if (EE < 0)
                EE = D20 * PI + EE;
            double M = EE - ES;
            double PZ = EC * Z / E / R - Math.Sqrt(A) * ES * ZD / E;
            double QZ = (ES * Z / E / R + Math.Sqrt(A) * (EC / E - E) * ZD) / Math.Sqrt(D10 - E * E);
            double RX = (Y * ZD - Z * YD) / Math.Sqrt(A * (D10 - E * E));
            double RY = (Z * XD - X * ZD) / Math.Sqrt(A * (D10 - E * E));
            double RZ = (X * YD - Y * XD) / Math.Sqrt(A * (D10 - E * E));
            double W = Math.Atan2(PZ, QZ);
            double L = Math.Atan2(RX,-RY);
            double I = Math.Acos(RZ);
            if (W < 0)
                W = D20* PI + W;
            if (L < 0)
                L = D20* PI + L;
            ZZ[0] = A * Ae;
            ZZ[1] = E;
            ZZ[2] = I;
            ZZ[3] = L;
            ZZ[4] = W;
            ZZ[5] = M;
            return ZZ;
        }

        public SateliteData calSpt(List<double> PVcur)
        {
            SateliteData sateliteData = new SateliteData();
            double We = 0.000072921150;
            double GE = 398600.4418e9;
            DateTime dateTime = new DateTime(y,m,d,h,mm,ss);
            DateTime date = new DateTime(2020, 1, 1, 0, 0, 0);
            DateTime Time = new DateTime(2006, 1, 1, 0, 0, 0);
            sateliteData.utcSecond = ModifiedJulianDate(dateTime) * 86400 - ModifiedJulianDate(date) * 86400;
            sateliteData.sptSecond = ModifiedJulianDate(dateTime) * 86400 - ModifiedJulianDate(Time) * 86400 + 4;
            PVcur[3] = PVcur[3] - PVcur[1] * We;
            PVcur[4] = PVcur[4] - PVcur[0] * We;
            List<double> dElementRef = PotoSig(PVcur);
            if (dElementRef[3] > Math.PI)
                dElementRef[3] = dElementRef[3] - 2 * Math.PI;
            if (dElementRef[4] > Math.PI)
                dElementRef[4] = dElementRef[4] - 2 * Math.PI;
            if (dElementRef[5] > Math.PI)
                dElementRef[5] = dElementRef[5] - 2 * Math.PI;
            double a_ref = 6500000;
            double SemiMajorAxis = dElementRef[0];
            double SemiMajorAxisOffset = SemiMajorAxis - a_ref;
            double SemiMajorAxisOffset1 = SemiMajorAxisOffset / 4.249 * Math.Pow(10, -3);
            double Eccentricity1 = dElementRef[1]/1.431*Math.Pow(10, -8);
            double Inclination1 = dElementRef[2] / Math.Pow(2, -27) / Math.PI;
            double RightAscensionOfAscendingNode1 = dElementRef[3] / Math.Pow(2, -27) / Math.PI;
            double ArgumentOfPerigee1 = dElementRef[4] / Math.Pow(2, -27) / Math.PI;
            double MeanAnomaly1 = dElementRef[5] / Math.Pow(2, -27) / Math.PI;
            sateliteData.delta_a = Math.Round(SemiMajorAxisOffset1);
            sateliteData.e = Math.Round(Eccentricity1);
            sateliteData.delta_i = Math.Round(Inclination1);
            sateliteData.OMEGA0 = Math.Round(RightAscensionOfAscendingNode1);
            sateliteData.omega = Math.Round(ArgumentOfPerigee1);
            sateliteData.M0 = Math.Round(MeanAnomaly1);
            return sateliteData;
        }

        public static double ModifiedJulianDate(DateTime date)
        {
            // 计算儒略日
            double julianDate = date.ToOADate() + 2415018.5;

            // 计算修正儒略日
            double modifiedJulianDate = julianDate - 2400000.5;

            return modifiedJulianDate;
        }
    }
}
